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Abstract. 

The Ensemble Kalman filter (EnKF) was introduced by Evensen in 1994 [1] as 
a novel method for data assimilation: state estimation for noisly observed time- 
dependent problems. Since that time it has had enormous impact in many application 
domains because of its robustness and ease of implementation, and numerical evidence 
of its accuracy. In this paper we demonstrate that the EnKF can be generalized to a 
wide class of inverse problems. In this context we show that the the EnKF estimate 
of the unknown function lies in a subspace A spanned by the initial ensemble. Hence 
the resulting error may be bounded above by the error found from best approximation 
in this subspace. We provide numerical experiments which compare the error incurred 
by EnKF for inverse problems with the error of the best approximation in A, and 
with variants on traditional least-squares approaches, restricted to the subspace A. 
In so doing we demonstrate that the EnKF method for inverse problems provides a 
derivative-free optimization method with comparable accuracy to that achieved by 
traditional least-squares approaches. Furthermore, we also demonstrate that accuracy 
is of the same order of magnitude as that achieved by best approximation. Three 
examples are used to demonstrate these assertions: inversion of a compact linear 
operator; inversion of piezometrix head to determine hydraulic conductivity in a Darcy 
model of groundwater flow; and inversion of Eulcrian velocity measurements at positive 
times to determine the initial condition in an imcompressible fluid. 



Submitted to: Inverse Problems 
1. Introduction 

Since its introduction in [1], the Ensemble Kalman filter (EnKF) has had enormous 
impact on applications of data assimilation to state estimation, and in particular 
in oceanography [2], reservoir modelling [3] and weather forecasting [4]; the books 
[5, 6, 7] give further details and references to applications in these fields. However, 
although demonstrably accurate, robust and simple to implement, the EnKF is not 
widely analyzed and its properties are poorly understood. The purpose of this paper 
is threefold: (i) to demonstrate that a novel non-standard perspective on the EnKF 
methodology converts it into a generic tool for solving inverse problems; (ii) to provide 
some basic analysis of the properties of this algorithm for inversion; (iii) to demonstrate 
numerically that the method can be effective on a wide range of applications. 
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In section 2 we describe the general inverse problem framework within which we 
work. We highlight the fact that the EnKF approach to this inverse problem provides 
an approximation to the solution in the subspace A spanned by the initial ensemble 
members, a fact observed for a specific sequential implementation of the EnKF in [8]. 
We then describe two classes of algorithms which we will use to evaluate the EnKF 
methodology. The first is best approximation of the truth in the subspace A. This is, 
of course, not an implementable method as the truth is not known; however it provides 
an important lower bound on the achievable error of the EnKF method and hence has 
a key conceptual role. The second class are variants on least squares, restricted to the 
subspace A. In section 3 we describe the EnKF method for the general inverse problem, 
which is in the form of a nonlinear iteration of an interacting ensemble. In Theorem 
3.1 we prove that the method gives rise to an approximation which lies in the subspace 
A spanned by the initial ensemble members and Corollary 3.3 uses this fact to give a 
lower bound on the achievable approximation error of the EnKF algorithm. 

Section 4 contains numerical experiments which illustrate the ideas in this paper 
on a linear inverse problem. The forward operator is a compact operator found from 
inverting the negative Laplacian plus identity. In section 5 we numerically study 
the groundwater flow inverse problem of determining hydraulic conductivity from 
piezometric head measurement in an elliptic Darcy flow model. Section 6 contains 
numerical results concerning the problem of determining the initial condition for the 
velocity in a Navier-Stokes model of an incompressible fluid; the observed data are 
pointwise (Eulerian) measurements of the velocity field. 

The numerical results in this paper all demonstrate that the EnKF method 
for inversion is a derivative-free regularized optimization technique which produces 
numerical results similar in accuracy to those found from least-squares based methods 
in the same subspace A. Furthermore, the three examples serve to illustrate the point 
that the method offers considerable flexibility through the choice of initial ensemble, 
and hence the subspace A in which it produces an approximation. In particular, for 
the linear and Darcy inverse problems, we make two choices of initial ensemble: (i) 
draws from a prior Gaussian measure and (ii) the Karhunen-Loeve basis functions of 
the centered Gaussian measure shifted by the mean. For the Navier-Stokes inverse 
problem the initial ensemble is also chosen to comprise randomly drawn functions on 
the attractor of the dynamical system. 

2. The Problem 

2.1. The Set- Up 

Let Q : X — > Y where X and Y are Hilbert spaces. In the following, we use (•, •) and 
|| • ||, etc. as the inner-product and norm on both X and Y, and it will be clear from 
the context which space is intended. Also, we let || • ||a denote \\A~^ ■ || for any positive 
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self-adjoint A. Consider the following equation: 
y = Q{u) +T]. 



(1) 



Here u G X is an unknown function, Q is the forward response operator mapping the 
unknown to the response space, r] G Y is a noise and y G Y the observed data. We 
assume that rj is a centered Gaussian with covariance T. The objective is to estimate 
u from y. We are particularly interested in situations where inversion of Q on Y is 
ill-posed. 

2.2. The Algorithms 

We will show that the EnKF method for inverse problems, defined in detail in the 
next section, provides an estimate of u which lies in the subspace A spanned by the 
initial ensemble used to seed the algorithm. Via an iterative interacting particle system, 
initiated at this ensemble, the method provides solutions of the inverse problem which lie 
in the subspace A and are formed of combinations of the initial ensemble with nonlinear 
weights reflecting the observed data. 

We wish to evaluate the EnKF method. Central to this task is the misfit functional 



which should be minimized in some sense. When inversion of Q is ill-posed, minimization 
of $ over the whole space X is not possible and regularization is required. For this paper 
three important regularization techniques [9] will play a role. The first is Tikhonov- 
Phillips regularization, whereby minimization of $ is replaced by minimization of the 
regularized functional 



where C is a compact, positive and self-adjoint operator on X and u and element of 
X. The second is regularization via minimization of $ over a compact set A G X. 
In this paper A will always be the linear span of a finite set of elements in X. The 
third is regularization via truncated iteration. The first, second and third methods 
can all be, and will be, used in conjunction. We note that minimization of I given by 
(3) has a statistical interpretation as the MAP estimator for Bayesian solution to the 
inverse problem with Gaussian prior N(u,C) on u [ )]. Noting this intepretation, we 
observe that we will use as choices for A both a set of draws from the Gaussian prior 
/io = N(u,C) and a subset of the the eigenf unctions of C - the Karhunen-Loeve basis 
- added to u. Finally we observe that, given the subspace A, and given the truth v) 
underlying the observed data y, we can compute the best approximation to the truth in 
A. Since the EnKF method presented here finds solutions confined to A it is natural to 
compare it with the best approximation within A. We now give more details of these 
least squares and best approximation algorithms used to benchmark the EnKF method. 
Furthermore, the following serves to establish a common notation used throughout our 
numerical studies. 




(2) 




(3) 
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Consider the space A = spa.n{ip^}j =1 comprised of the initial ensemble used for 
the EnKF. Let fiQ be some prior measure. The space A will be chosen based on the 
prior /i . When /i is Gaussian N(u,C) let (Xj,4>j) denote eigenvalue/eigenvector pairs, 
in descending order by eigenvalue - this is the Karhunen-Loeve basis. However we do 
not restrict ourselves to Gaussian priors: for the Navier-Stokes example the prior will 
be the measure supported on the attractor. We will use: (a) ip^) ~ ^ i.i.d. and all 
algorithms using this choice of A will bear the subscript R for random; for Gaussian 
priors, we will additionally use (b) ij)^> — u + a/A^-, j < J, and all algorithms using 
this method will bear the subscript KL for Karhunen-Loeve . 

We will compute EnKF^ and EnKFxL approximations to the inverse problems of 
interest. For comparison we will compute the best approximation 

u B a = axgmin ue4 ||u - v?\\l. (4) 

leading to BAr and BAkl- Then, we will compute the projected basis least squares 
solution 

uls = argmin ue J \y - G(u) 1 1£ (5) 

or generalizations to include Tikhonov-Phillip regularization [11], based on fiQ = 
N(u,C), and truncated Newton-CG iterative methods [12]. In all cases we will denote 
the resulting approximations denoted by LSr and LSkl- 

2.3. Linear Problems 

It is instructive to consider the case where Q{u) = Gu for some linear operator 
G : X — > Y as this will enable us to makes links with standard regularized least squares 
problems, the Kalman filter and hence the ensemble Kalman filter. If \\y — G ■ ||p is 
continuous on D(C) then the Tikhonov-Phillips regularized solution from (3) is given 



by mtp = argmin ug £> (c) /(w), where 

I{u) = \\y — Gu\\y + \\u — u\\q. (6) 

Furthermore this minimization has a unique solution which and may be computed 
explicitly to be 

u TP =u + CG*(GCG* + T)-\y-Gu). (7) 

To see this we apply the Sherman- Morrison- Woodbury formula [13] to the normal 
equations resulting from minimization of /, as follows: 

u TP = (GT^G + C- 1 ) ^ (G*T^ 1 y + C^u) (8) 

= u+{G*T- l G + C- l y 1 G*T-\y-Gu) (9) 

= u+[C- CG*(T + GCG*Y l GC] G*T-\y - Gu) (10) 

= u + [CG* - CG*(T + GCG*Y l GCG*] T-\y - Gu) (11) 

= u + [CG* - CG* (I - (r + GCG*)-^)] r- x (y - Gu) (12) 

= u + CG*(T + GCG*)- 1 {y-Gu), (13) 
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which we see immediately is exactly equation (7). 

It is well-known that the Tikhonov regularized minimizer is the mean of the 
conditional random variable u\y [14, 15], where the joint random variable is prescribed 
by the identities 

u~N(u,C) (14) 
y = Q u + r]i v „ jV(0,r) ±u (15) 

and where _L denotes independence. Furthermore, we can connect this conditional 
random variable with the Kalman filter [16] as we now show. Let Z = X xY and define 
z G Z and the linear operator A : Z — > Z by 

Az = 

Define the operator H : Z -»■ Y by H = (0, 1), and H x : Z -»■ X by H 1 - = (I, 0). Now 
consider the partially observed linear Gaussian dynamical system given by 

Vn+1 = Hz n+X + r] n +i. 

Here it is assumed that {r] n } n£ z+ is an i.i.d sequence with 771 ~ N(0,T). The Kalman 
filter computes the Gaussian random variable z n \{y m }'!^ n=l , which is characterized by its 
mean and covariance for each n, and can be implemented sequentially as an iteration 
in n. If uq ~ N(u, C), po = Guo and yi — y then, by construction, the first step of the 
Kalman filter Z\\yi, projected into the u coordinate by H x , is identical to the random 
variable u\y defined previously. In particular the mean of one step of this Kalman filter 
is equal to utp- Note, further, that this link between the linear inverse problem suggests 
the intruiguing possibility of iterating the Kalman filter to get information about the 
inverse problem. 

The EnKF proceeds by moving an ensemble of particles according to the mean 
Kalman update, and using this ensemble to approximate the covariance information 
needed to update the mean. In the limit of an infinite ensemble, EnKF converges to 
the Kalman Filter for linear problems. As we will show in section 3, when applied 
to a general class of nonlinear inverse problems, including the linear example under 
consideration here, the EnKF produces an approximation in the linear span of the initial 
ensemble. As such there is a natural connection to Krylov methods [17] and minimum 
residual methods such as MINRES [18] and GMRES [19] for appropriate choices of this 
initial ensemble. 

Notice that we have made two interesting observations, in the context of this steady 
linear inverse problem: (i) we have connected it to a dynamic iterative estimation 
scheme; (ii) we have remarked that covariance information needed to implement the 
iteration can be approximated by an ensemble. These ideas can be generalized to 
nonlinear problems, and this is the subject of the next section. 
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3. EnKF Algorithm 

The idea of the Kalman filter may be extended to nonlinear problems via the Extended 
Kalman Filter (ExKF) [: )] which proceeds by (i) approximating the forward model by its 
linearization Q{u) ~ Q{u) +DQ{u){u — u) and (ii) assuming that the desired distribution 
is well-approximated by a Gaussian u ~ N(u, C). It is then possible to apply the Kalman 
Filter, using this approximation. The EnKF [1] makes the further approximation of 
working with an ensemble and using this to estimate covariance information. For 
sequential problems numerical evidence demonstrates that the EnKF can produce 
reasonable state estimation even outside the regimes where the approximations (i) and 
(ii) are valid - see [21]. In this section we show how the EnKF can be used to construct 
an algorithm for inverse problems, by introducing an iterative scheme, generalizing 
subsection 2.3. 

Define Z = X x Y , z & Z, H : Z — > Y and H 1 - : Z — > X as in the previous 
subsection. Then define 5 : Z — > Z by 



We then estimate u by a variant of the EnKF algorithm applied to a partially observed 
dynamical system of the form 



Here it is assumed that {r] n } n& z+ is an i.i.d. sequence with r/i ~ N(0,T). The EnKF 
estimates z n given {y m }^ =1 . Since H 1 - z = u we will apply the projector H 1 - to our 
estimate of the state of the system at the n th iterate of the algorithm in order to 
determine the unknown function u. 

Recall that y is the given data. We will construct articifical data {y n }n&+ from 
this single observation y e Y in two ways, both described below. The algorithm works 
with an ensemble of vectors {zn }j=i which represent an ensemble of estimates of the 
state of the system at the n th step of the algorithm, given the data {y m }m=i- We will 
take the mean of this ensemble as an estimate of the state. We assume that z ( j) = ip^ 
for I < j < J. The set A = span{-0(^}^ =1 will play a key role in the analysis of the 



The algorithm breaks into two parts, a prediction step and an analysis step, and 
we desribe these separately. We then discuss properties of the resulting algorithm. The 
prediction step maps the current ensemble of particles into the data space, and thus 
introduces information about the forward model. The analysis step makes comparisons 
of the mapped ensemble, in the data space, with the data, or with noisily perturbed 
versions of the data; it is at this stage that the ensemble is modified in an attempt to 
better match the data. 




y n+ i = Hz n+ i + r] n+1 . 



z., 



(16) 
(17) 



algorithm. 
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3.1. Prediction Step 

In this step we propagate the ensemble of particles forward under (16) giving 

tli = (18) 
From this ensemble we define a sample mean and covariance as follows: 



4 + i=7ES ( 19 ) 



Cn+1 = -J 2_j Z n+l( Z n+l) T ~~ Z n+l Z n+V (20) 

i=i 

All the vectors and operators involved have block structure inherited from the structure 
of the space Z = X x Y. For example we have 

J01 _ ( \ _ ( »" W \ M) _ f ^ 



We also have 



The vectors u n and p n are given by 



j=l 3=1 



Pn = 7E^' ) = 7E^ 1 ) (22) 



The blocks within C n are given by 
j 



i=i 

C7= 1^^) ( ^) ) T _- n -T (24) 



j 



3.2. Analysis Step 



\Y.f^lg*) T -pJ>Z,. (25) 



The analysis step combines the ensemble {^+i}/ = i with data {y„+i}/=i, using ideas 
stemming from the Kalman filter, but with the predicted mean and covariance estimated 
from the samples. For the perturbed observation variant of EnKF the noisly perturbed 
data is constructed as follows: 

V$i = V + V$.v (26) 
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Here the rftt+i are an i.i.d collections of vectors indexed by (j,n) with rj^ ~ N(0,T). 
We will also study, for comparison, an unperturbed algorithm where we simply take 

y% = y- (27) 

We do not include the results for this they are sufficiently similar to the results 

for the case of perturbed observations. 

To construct the analysis step we define the operators S n and K n by 

S n = HC n H T + T, (28) 
K n = C n H T S-\ (29) 



where 



H 1 = 




The operator K n+1 is the Kalman gain and the particles are updated according to the 
standard Kalman update formula which blends model and data: 



z 



n+l 



I^li + KnMU-Hz^l,) (30) 



= (/ - K n+1 H)^ +1 + K n+l y%. (31) 
3.3. Properties of the Algorithm 

Recall the space A formed from the linear span of the initial ensemble. All subsequent 
EnKF vectors remain within this space. 

Theorem 3.1. For every (n,j) e N x {1, • • • , J} we have Zn^ G A. 

Proof. This is a straightforward induction, using the properties of the update formulae. 
Clearly the statement is true for n — 0. Assume that it is true for n. The operators S n 
and K n have particular structure inherited from the form of H. In particular we have 



s n = cr + r 



and 

Note that 

and 

where 



crtcr + r)- 1 



cS p = 7EpS?" ) (p? ) )' j 



1 3 



= 1 
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Recall that, from the structure of the map S, we have 

u n+l — a n 1 Pn+l — =A u n ) 

and from the definition of H we get that 

o c^cw + r)- 1 
o epic** + T)- 1 



Using these facts in (30) we deduce that the update equations are 

= »%> + CZxV&x + r)- 1 (*& - ( 32 ) 
Si = G(rfP) + cZAcZi + r)- 1 - • (33) 



Pr. 

If we define 



then the update formula (32) for the unknown u may be written as 



u% = ug> + j ^(pit, (34) 



fc=i 

k=l 

In view of the inductive hypothesis at step n, this demonstrates that u^ Srl G A for all 
j G {1, • • • , J} and the proof is complete. □ 

Remark 3.2. The proof demonstrates that the solution at step n is simply a linear 
combination of the original samples; however the coefficients in the linear combination 
depend nonlinearly on the process. Note also that the update formula (33) for the system 
response may be written as 

pSi = + j OpSi (36) 

fc=i 
j 



= Qi^) + j £<pSi, d%)g(u^). (37) 
fc=l 

Let denote the true unknown function which underlies the data so that 

y = Q{u ] ) + 

for some noise G Y. We have the following lower bound for the accuracy of the 
estimates produced by the EnKf algorithm for inverse problems: 

Corollary 3.3. The error between the estimate of u at time n and the truth v) satisfies 



u n -u 



t|| > inf \\v — u*\ 



veA 
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4. Elliptic Equation 

As a simple pedagogical example, we consider the ill-posed inverse problem of recovering 
the right-hand side of an elliptic equation in one spatial dimension given noisy 
observation of the solution. This explicitly solvable linear model will allow us to 
elucidate the performance of EnKF as an iterative regularisation method. The results 
show that EnKF performs comparably to Tikhonov-Phillips regularized least squares 
for this problem. The lower bound BA is always significantly better, but of course 
unimplementable in practice, since the truth is unknown. 

4-1. Setting 

Consider the one dimensional elliptic equation 



Thus G = A' 1 where A = + 1) and D(A) = H 2 (I) n H%(I) with I = (0, tt). We 



are interested in the inverse problem of recovering u from noisy observations of p: 



For simplicity we assume that the noise is white: 77 ~ N(0, 7 2 /). We consider Tikhonov- 
Phillips regular izat ion of the form where C = /3(A — The problem may be 
solved explicitly in the Fourier sine basis, and the coefficients of the Tikhonov-Phillips 
regularized least squares solution, in this basis, Uk may be expressed in terms of the 
coefficients of the data in the same basis, t/k'- 



This demonstrates explicitly the regularization which is present for wavenumbers k such 
that k 6 > 0(P , y~ 2 ). We now present some detailed numerical experiments which will 
place the EnKF algorithm in the context of an iterative regularization scheme. This 
will provide a roadmap for understanding the nonlinear inverse problem applications we 
explore in subsequent sections. 

4-2. Numerical Results 

Throughout this subsection we choose /3 = 10 and 7 = 0.01. We choose a truth 
v) ~ N(0,C) and simulate data from the model: y = Gu^ + rf, where ~ A/"(0,r). 
Recall from Section 2 that the space A = span{ip^}j =1 will be chosen either based on 
draws from N(0,C) (with subscript R) or from the Karhunen-Loeve basis for N(0,C) 
(with subscript KL.) 



d 2 p 

dx 2 ^ 



u(0) = u(tt) = 0. 



y = p + v 

= A~ x u + 77. 




k — 1, 00 



(38) 
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We evaluate the lower bounds on the error resulting from best approximation in 
the random and Karhunen-Loeve bases respectively: BAr and BAkl- We also compute 
errors in the projected basis regularized least squares solutions LSr and LSkl with 
additional Tikhonov-Phillips regularisation to compensate for the large number of basis 
vectors. And finally, we compute the errors resulting from the EnKF means, namely 
EhKFr and EuKFkl, as detailed in the previous section. Notice that for the EnKF 
methods we may take just one-step, or we may iterate over many steps; experiments 
describing both are reported. 

In Figure 1 (top) we show three examples of the error of the estimators with respect 
to the truth over 100 different ensembles A of random draws from the prior: one step of 
EnKFft, Tikhonov-Phillips regularized LSr (equivalent to KF) and BAr. The methods 
clearly indicate that the EnKF method is comparable to the LS method in terms of 
accuracy, but involves no matrix inversions. Both EnKF and LS are less accurate than 
BA, of course, but produce errors of similar order of magnitude. 

In Figure 1 (bottom) the results for random draws, with suffix R, are averaged 
over many initial choices of ensembles, and compared with results using the Karhunen- 
Loeve basis, with suffix KL. The errors with suffix R are therefore averaged over 
randomly drawn initial bases and horizontal lines indicate these quantities in the plot. 
EnKF based methods are compared with BA and LS type methods. These results again 
show that EnKF and LS are similar in terms of accuracy, and produce errors of similar 
order of magnitude to BA. Furthermore they show that, when iterated, the error of 
EnKFxL decreases for some time before reaching its minimum, while EnKF^ reaches its 
minimum after a small number of iterations (often only one), and then increases (note 
that the averaging of the error hides the actual value at which the error starts increasing 
for any given choice of initial ensemble). 

5. Groundwater flow 

Next, we will investigate an inverse problem arising in reservoir modeling. Typically, 
the calibration of subsurface flow models consist of estimating geologic properties whose 
model predictions "best fit" the flow measurements [22]. In other words, the calibration 
can be posed as the inverse problem of finding the geologic properties given subsurface 
flow data. In particular, herein we will consider the estimation of the conductivity of 
an aquifer form hydraulic head measurements. Similarly to the previous section, we 
find comparable performance of EnKF and LS, with the latter here regularized in a 
Newton-CG fashion. Again, BA is included for comparison. 

5.1. Setting 

We consider groundwater flow in a two-dimensional confined aquifer whose physical 
domain is O = [0, 6] X [0, 6]. The hydraulic conductivity is denoted by K. The flow in 
the aquifer, is described in terms of the piezometric head h(x) (x G Q) which, in the 
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Figure 1. Top: Comparison over an ensemble of random ensembles A of the relative errors 
of one iteration of EuKFr versus (Tikhonov regularized) least squares (LSr,) and the best 
approximation (BAr,). Bottom: Comparison of one trajectory of EnKFp, and EuKFkl over 
iteration n, as well as EuKFr, BAr and LSr averaged over ensembles, and BAkl and LSkl- 



steady-state is governed by the following equation [23] 
- V • e u Vh = f in fi 



(39) 
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where u = \ogK and / is defined by 

{0 if < x 2 < 4, 
137 if 4 < x 2 < 5, (40) 
274 if 5 < x 2 < 6. 

We consider the following boundary conditions 

fc(x,0) = 100, ^-(6,y) = 0, -e^(0,y)= 500, ™{x,6) = 0, (41) 

For the physical interpretation of the source term and boundary conditions (40)-(41), 
we refer the reader to [ ] where a similar model was used as a benchmark for inverse 
modeling in groundwater flow. A similar model was also studied in [25, 26] also in 
the context of parameter identification. We will be interested in the inverse problem 
of recovering the hydraulic conductivity, or more precisely its logarithm it, from noisy 
pointwise measurements of the piezometric head h. This is a model for the situation in 
groundwater applications where observations of head are used to infer the conductivity 
of the aquifer. 



5.2. Numerical Results 

We now explain in detail the numerical results presented in Figures 2 and 3. We let 
Q(u) = {p{xk)}keK where IK is some finite set of points in Q with cardinality K. In 
particular, K is given by the configuration of N = 100 observation wells displayed in 
Figure 3 (Top right). We introduce the prior Gaussian fi ~ J\f(u,C), let ~ fio, and 
simulate data y = G(vJ) + rf, where ry ~ Af(0,T). As in the previous section, we let 
C = (3(—A)~ a and T = 7 2 J, and here we choose a = 1.3, = 0.5, u = 4 and 7 = 7 fixed. 
We reiterate from Section 2 that the space A = span{^^)}^ =1 will be chosen based on 
either draws from the prior fiQ, with subscript R, or on the Karhunen-Loeve basis, with 
subscript KL. 

The forward model (39)-(41) is discretized with cell-centered finite differences [27]. 
For the approximation of the LS problem, we implemented the Newton-CG method 
of [12]. We conduct experiments analogous to those of the previous section and the 
results are shown in Figure 2 and Figure 3. The results in Figure 2 are very similar to 
those shown in the previous section for the linear elliptiuc problem, demonstrating the 
robustness of the observations made in that section for the solution of inverse problems 
in general, using the EnKF methodology herein. 



6. Navier-Stokes Equation 

In this section, we consider an inverse problem in fluid dynamics, which is relevant 
to data assimilation applications in oceanography and meteorology. In particular, we 
examine the problem of recovering the initial condition of the Navier-Stokes Equation 
(NSE), given noisy pointwise observations of the velocity field at later times. We will 
investigate a regime in which the combination of viscosity, time-interval, and truncation 
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Figure 2. Top: Comparison over an ensemble of random ensembles A of the relative errors 
of one iteration of EnKFR versus least squares (LSr) and the best approximation (BAr). 
Bottom: Comparison of one trajectory of EnKFp and EnKFxL over iteration n, as well as 
EnKFR, BAr and LSr averaged over ensembles, and BAkl and LSrl- 



of the forward model is such that the inverse problem is only "mildly ill-posed" , in the 
sense that it does not exhibit the ill-posedness of the function space limit. 
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Figure 3. Estimated log K. Top left: BAr, average. Top middle left: LSr, average. 
Top middle right: first iteration of EnKFR, average. Top right: measurement well locations. 
Bottom left: BAkl- Bottom middle left: LSkl- Bottom middle right: EnKFxL- Bottom 
right: the truth . 

In practice, a regularization term is included anyway, but with the purpose of 
incorporating background, or climatological, information into the model in the form of 
a prior. 

As in the previous two sections we will see that EnKF and LS type methods 
perform comparably, in both the case of random and Karhunen-Loeve based initial 
draws. Furthermore the lower bound produced by BA type methods is of similar order of 
magnitude to the EnKF-based methods although the actual error is, of course, smaller. 
However, we will also see that the behavior of the iterated EnKFR, is quite different from 
what we observed in the previous sections. We conjecture that this is because of the 
mildly ill-posed nature of this problem. Indeed we have have repeated the results of this 
section at higher viscosity, where linear damping in the forward model induces greater 
ill-posedness, and confirmed that we recover results for the iterated EnKFR which are 
similar to those in the previous sections. 

6.1. Setting 

We consider the 2D Navier-Stokes equation on the torus T 2 := [—1, 1) x [—1, 1) with 
periodic boundary conditions: 

d t v - uAv + v-Vv + Vp = f for all (x, t) e T 2 x (0, oo), 
V • v = for all (x, i)eT 2 x (0, oo), 

v = u for all (x,t) G T 2 x {0}. 

Here v. T 2 x (0, oo) — > M 2 is a time-dependent vector field representing the velocity, 
p: T 2 x (0, oo) — > R is a time-dependent scalar field representing the pressure, /: T 2 — > IR 2 
is a vector field representing the forcing (which we assume to be time-independent 
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for simplicity), and v is the viscosity. We are interested in the inverse problem of 
determining the initial velocity field u from pointwise measurements of the velocity field 
at later times. This is a model for the situation in weather forecasting where observations 
of the atmosphere are used to improve the initial condition used for forecasting. 

6.2. Numerical Results 

We let tj = jh, for j = 0, . . . , J, and Q{u) = {v(xk, tj)}(j,k)ew where K' — {0, • • • , J} x K 
and IK is some finite set of points in Q with cardinality K. In particular, we take K to 
be the set of grid points in physical space implied by the underlying spectral truncation 
used in the numerical integration (details below). As discussed in Section 2, we introduce 
a prior /i , let ~ fi , and simulate data y = Qiv)) + , where ~ A^(0,r). As in 
the previous section, we let T = 7 2 /, so our tunable parameter is 7 which we fix at 
7 = 0.01 for these experiments. The prior /i is defined to be the empirical measure 
supported on the attractor, i.e. it is defined by samples of a trajectory of the forward 
model after convergence to statistical equilibrium. We use 10 4 time-steps to construct 
this empirical measure. Once again the space A = span{?/'^)}^ =1 is chosen based on J 
samples from the (now non-Gaussian) prior no, or on a Karhunen-Loeve expansion based 
on the empirical mean u and covariance C from the long forward simulation giving rise 
to /io- It is important to note that the truth is included in the long trajectory used to 
construct fiQ. Therefore, some of the initial ensembles drawn from ^ en d up containing 
a snapshot very close to the truth; such results are overly optimistic as this is not the 
typical situation. It is interesting that each method has comparably overly optimistic 
results for such ensembles. For (ii), we use the the empirical mean u and covariance C 
over a long forward simulation. 

The forcing in / is taken to be / = V~ L, ?/S where if) = cos(irk ■ x) and V -1 = JV 
with J the canonical skew-symmetric matrix, and k = (5,5). The method used to 
approximate the forward model is a modification of a fourth-order Runge-Kutta method, 
ETD4RK [ 8], in which the Stokes semi-group is computed exactly by working in 
the incompressible Fourier basis {4>k(x)}k<=z 2 \{o}, an d Duhamel's principle (variation 
of constants formula) is used to incorporate the nonlinear term. We use a time-step of 
dt = 0.005. Spatially, a Galerkin spectral method [ )] is used, in the same basis, and the 
convolutions arising from products in the nonlinear term are computed via FFTs. We 
use a double-sized domain in each dimension, padded with zeros, resulting in 64 2 grid- 
point FFTs, and only half the modes in each direction are retained when transforming 
back into spectral space again. This prevents aliasing, which is avoided as long as more 
than one third of the domain in which FFTs are computed consists of such padding 
with zeros. The dimension of the attractor is determined by the viscosity parameter v. 
For the particular forcing used there is an explicit steady state for all v > and for 
v > 0.035 this solution is stable (see [30], Chapter 2 for details). As v decreases the flow 
becomes increasingly complex and we focus subsequent studies of the inverse problem 
on the mildly chaotic regime which arises for v = 0.01. Regarding observations, we 
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let h — 4 x dt — 0.02 and take J = 10, so that T = 0.2. The observations are made 
at all numerically resolved, and hence observable, wavenumbers in the system; hence 
K = 32 2 , because of the padding to avoid aliasing. 

The numerical results resulting from these experiments are displayed in in Figures 
4 and 5. The results are very similar to those of the previous two sections, qualitatively: 
EnKF and LS type methods perform comparably, in both the case of random and 
Karhunen-Loeve based initial draws; furthermore the lower bound produced by BA 
type methods is of similar order of magnitude to the EnKF-based methods although 
the actual error is, of course, smaller. However, we also see that the behavior of the 
iterated EnKFp is quite different from what we observed in the previous sections since it 
decreases montonically. We conjecture that this is because of the mildly ill-posed nature 
of this problem. We also note that Tikhonov of truncated iterative regularisation is not 
necessary for LS for this example, because of the truncation of wavenumbers, the low 
viscosity, and the short time-interval. We have have repeated the results of this section 
at the higher viscosity v = 0.1, where linear damping in the forward model induces 
greater ill-posedness, and confirmed that we recover results for the iterated EnKF^ 
which are similar to those in the previous sections, in that the error displays a minimum 
after one or two steps Also, for this higher value of u, Tikhonov regularization is used to 
ensure that LS delivers reasonable results, because of the greater ill-posedness. Finally, 
we have checked that the behavior of the error in the iterated EnKF^ is repeatable 
for Gaussian prior /z , and hence is not a result of non-Gaussian ensembles used in the 
figures. 

7. Conclusions 

We have illustrated the use of EnKF as a derivative-free optimization tool for inverse 
problems, showing that the method computes a nonlinear approximation in the linear 
span of the initial ensemble. We have also demonstrated comparable accuracy to least 
squares based methods and shown that, furthermore, the accuracy is of the same order of 
magnitude as best approximation within the linear span of the initial ensemble. Further 
study of the EnKF methodology for inverse problems, and in particular its accuracy with 
respect to choice of initial ensemble, would be of interest. Furthermore, in this paper 
we have concentrated purely on the accuracy of state estimation using EnKF. Study 
of its accuracy in terms of uncertainty quantification will yield further insight. Finally, 
although we have studied a time-dependent example (the Navier-Stokes equation) we 
did not use a methodology which exploited the sequential aquisition of the data - we 
concatenated all the data in space and time. In future work it is of interest to study 
ideas similar to those herein, but exploiting sequential structure. 
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Figure 4. Top: Comparison over an ensemble of random ensembles A of the relative errors 
of one iteration of EnKF^ versus least squares (LSr) and the best approximation (BAr). 
Bottom: Comparison of one trajectory of EnKFfj, and EnKFxL over iteration n, as well as 
EuKFr, BAr and LSr averaged over ensembles, and BAkl and LSkl- 
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